package scu.maqiang.elasticity;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T32;
import scu.maqiang.mesh.Mesh2T3;
import scu.maqiang.mesh.T3Type;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.DirectSSolver;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.SRMatrix;

public class Elasticity2DT3Case1 {
    public static void main(String[] args) {
        Mesh2T3 mesh = new Mesh2T3().square2D(40, 4, T3Type.Left);
        mesh.scale(10.0, 1.0);

        FES2T32 fsElastic = new FES2T32(mesh);
        int dof = fsElastic.GetNdof();
        SRMatrix K = new SRMatrix(dof);
        double[] RHS = new double[dof];
        double[] ENu = {2.0e6, 0.3};
        fsElastic.assembleStiff(ENu, BVPType.PLANESTRAIN, K);
        System.out.println(K);
        double[] force = {0.0, -1.0};
        fsElastic.assembleSource(force, BVPType.COMMON, RHS);
        //System.out.println(MVO.toString(RHS));
        fsElastic.applyBC_MBN(K, RHS, Direct.All, 0.0, new int[]{4});
        DirectSSolver solver = new DirectSSolver(K);
        double[] sol = solver.CholeskySolve(RHS, 1);
        double[][] disp = new double[2][mesh.getNv()];
        fsElastic.extractComponent(sol, disp);

    }
}
